program SimObsComp

    use iasing_osse, only:  &
        JPI4,               &
        JPR4,               &
        JPR8,               &
        NIASI,              &
        NLM60,              &
        NIANG,              &
        r4a_ids,            &
        r60_dat,            &
        GAS_ID_CO2,         &
        GAS_UNIT_SPECCONC

    use netcdf

    implicit none
    integer :: jo
    ! 
    integer(JPI4)       :: dbid
    type(r4a_ids)       :: db_ids
    type(r60_dat)       :: db_dat
    character(len=128)  :: dbfile

    logical             :: lbrt, lex

    integer(JPI4)       :: natm, nocl, nlev
    integer(JPI4), allocatable  :: ocl(:)

    real(JPR4), dimension(NIASI)    :: wia
    real(JPR4), allocatable         :: spc(:,:), obt(:,:), co2(:,:), q(:,:)
    real(JPR4), allocatable         :: s1(:), s2(:)
    integer(JPI4), allocatable      :: sn(:)
    real(JPR4), dimension(NLM60)    :: lev, co21, co22, co2n
    integer :: jl
    real(JPR8) :: dpaux

    include "checknc.interface"
    include "r60_initializer.interface"
    include "r60_reader.interface"
    include "rad2brt_f.interface"
    include "rttov_q2v.interface"

    data lbrt /.false./

    !dbfile = "/home/javier/CurrentWork/00_Data/simul_rt.nc"
    
    call get_command_argument(1,dbfile)
    inquire(file=dbfile, exist=lex)
    if (.not. lex) then
        print*, "File '", trim(dbfile), "' does not exists. Aborting execution..."
        stop -1
    endif

    call checknc(nf90_open(dbfile, NF90_NOWRITE, dbid))
    db_ids = r60_initializer(dbid)

    ! Inquire
    call checknc(nf90_inquire_dimension(dbid, db_ids%atm_did, len=natm))
    call checknc(nf90_inquire_dimension(dbid, db_ids%ocl_did, len=nocl))
    call checknc(nf90_inquire_dimension(dbid, db_ids%lev_did, len=nlev))
    !print*, natm, nocl, nlev
    !print*, NIASI*nocl*4./1024**2

    ! Allocate arrays
    allocate(spc(NIASI,natm))
    allocate(obt(nocl,natm))
    allocate(co2(NLM60,natm))
    allocate(q(NLM60,natm))
    allocate(ocl(nocl))

    ! Allocation of statistical arrays
    allocate(s1(nocl))
    allocate(s2(nocl))
    allocate(sn(nocl))

    ! Load data
    call checknc(nf90_get_var(db_ids%dim%gid, db_ids%dim%wia, wia))
    call checknc(nf90_get_var(db_ids%dim%gid, db_ids%dim%ocl, ocl))
    call checknc(nf90_get_var(db_ids%rbt%gid, db_ids%rbt%ria, spc))
    call checknc(nf90_get_var(db_ids%obs%gid, db_ids%obs%obt, obt))

    ! Change units if required
    if (spc(1,1) < 100. .and. spc(1,1)>0.) then
        lbrt = .true.
        do jo = 1, natm
            spc(:,jo) = rad2brt_f(NIASI, wia, spc(:,jo))
        enddo
    endif
    if (lbrt) then 
		print*, "Units Changed"
		!stop
    endif


	sn = count(spc>0., dim=2)
    s1 = sum(spc(ocl,:) - obt, dim=2, mask=spc>0) / sn
    s2 = sqrt((sum((spc(ocl,:) - obt)**2,dim=2, mask=spc>0) &
        - s1**2*sn)/(sn-1))

    !print '(10F12.6)', spc(ocl, 1) - obt(:,1)
    
    do jo = 1, nocl
        write(*,'(I6,F10.3,2F10.2,I10)') ocl(jo), wia(ocl(jo)), s1(jo), s2(jo), sn(jo)
    enddo

    ! Arrays deallocation 
    deallocate(obt)
    deallocate(ocl)
    deallocate(spc)

    deallocate(s1)
    deallocate(s2)
    deallocate(sn)

    ! Close file
    call checknc(nf90_close(dbid))

end program
